Solution of large scale nuclear structure problems by wave function factorization 
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Low-lying shell model states may be approximated accurately by a sum over products of proton 
and neutron states. The optimal factors are determined by a variational principle and result from 
the solution of rather low-dimensional eigenvalue problems. Application of this method to sd-shell 
nuclei, p/-shell nuclei, and to no-core shell model problems shows that very accurate approximations 
to the exact solutions may be obtained. Their energies, quantum numbers and overlaps with exact 
eigenstates converge exponentially fast as the number of retained factors is increased. 
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I. INTRODUCTION 

Realistic nuclear structure models are difficult to solve 
due to the complexity of the nucleon-micleon interaction 
and the sheer size of the model spaces. Exact diagonal- 
izations are now possible for of -shell nuclei P, 0, iSj and 
for sufficiently light systems 0,13, and Quantum Monte 
Carlo calculations 0, Q have solved light nuclei up to 
about mass A = 12. For cases where an exact solution is 
not feasible, various approximations are employed. We 
mention stochastic methods like shell-model Monte Carlo 

and Monte Carlo shell-model > and recent ap- 
plications of coupled cluster expansions 0, . 

In recent years, several truncation method for shell 
model diagonalizations have been developed. These 
methods aim at a reduction of the enormous dimension- 
ality of shell model Hilbert spaces while maintaining a 
high accuracy in the computed observables. Based on 
arguments of statistical spectroscopy, Horoi and cowork- 
ers [l3,14,,15,, 16] developed the exponential convergence 
method. Mizusaki and Imada devised an extrap- 

olation method and applied it to a configuration trunca- 
tion. Both techniques use single-particle basis states and 
provide a method to extrapolate the results of truncated 
calculations to the full Hilbert space. Other approaches 
use correlated basis states to obtain a rapid convergence. 
Andreozzi et al. l20l |. for instance, construct a ba- 
sis from products of correlated proton states and corre- 
lated neutron states. Gueorguiev et al. use a 
mixed-mode shell model of single-particle and collective 
configurations, while Vargas et al. |23l | use a truncation 
based on coupled SU{'i) irreps to describe the interplay 
and competition of collective and single-particle degrees 
of freedom. 

Though the selection of the relevant states is physi- 
cally well motivated for all these truncation schemes, it 
does not directly follow from a variational principle. This 
is different for the density matrix renormalization group 
(DMRG) l24l and the very recently proposed factoriza- 
tion method [23 • The DMRG uses a sophisticated renor- 
malization and truncation scheme that includes the most 
important states and correlations. Dukelsky et al. p^l2^ 



and Dimitrova et al. |28j applied this method to nuclear 
structure problems. The ground-state factorization is 
based on a related truncation. At a given truncation, the 
optimal states are determined from a variational princi- 
ple. These last two methods also allow for an extrapola- 
tion to full Hilbert spaces as the results tend to converge 
exponentially. In this article we present a detailed de- 
scription of the factorization method and discuss several 
applications. A summary of some of the main results has 
been presented in an earlier paper psj . 

This article is organized as follows. In Section ^ we 
give a derivation of the main theoretical results and 
present details of the numerical implementation. Sec- 
tion lllll presents numerical results for sd-shell nuclei, pf- 
shell nuclei and for no-core shell-model problems. The 
convergence of the factorization method and a com- 
parison with other truncation methods is presented in 
Sect. II Vl and we conclude with Section Ivl 



II. THEORETICAL BACKGROUND 

A. Motivation 

Shell model basis states are products of proton Slater 
determinants {|7rQ),Q; = l,...,dp} and neutron Slater 
determinants {\va),a = l,...,fiAr}. Here, dp and 
denote the dimension of the proton space and the neutron 
space, respectively. The shell-model ground-state may be 
expanded as 



(1) 



This expansion is not unique since the amplitudes ^ ap 
depend on the choice of basis states within the two sub- 
sets. There is, however, a preferred basis in which the 
amplitudes '^ap are "diagonal" . This basis is formally 
obtained from a singular value decomposition of the rect- 
angular amplitude matrix '5 of dimension dp x d^ as 
= USV^. Here U {V) is a dp x dp {dn x dN) di- 
mensional matrix with orthonormalized columns. S* is a 
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FIG. 1: Singular values for ground states of ^°Ne {dp = 
66), 2^Ne {dp = 66), ^''Mg {dp = 495), ^**Mg {dp = 495), and 
^*Si {dp = 924). 

rectangular matrix of dimension dp x ^a? with elements 
S'ij = for i ^ j and the "diagonal" elements Si^i = Si 
are the non-negative singular values. Performing the sin- 
gular value decomposition yields 

min [dp .d^ ) 

l*)= E ^APo)\^o)- (2) 

Here Sj denote the singular values while the proton-states 
\pj) and the neutron-states \nj) are orthonormal sets of 
left and right singular vectors, respectively. In general, 
these states are superpositions of many Slater determi- 
nants and exhibit strong correlations. The non-negative 
singular values si > S2 > sz ■ ■ ■ fulfill ^"^ — ^ ^'^^ ^o 
wave function normalization. 

It is interesting to compute the singular value decom- 
position for ground states of realistic nuclear many-body 
Hamiltonians. To this purpose we perform a numeri- 
cal singular value decomposition of the amplitude matrix 
in Eq. (QJ using the Lapack routines [l^. Figure^] 
shows the squares of the singular values for the ground 
states of sd-shell nuclei ^°Ne, ^^Ne, ^''Mg and ^^Si (from 
the USD interaction [33) plotted versus their normalized 
index j/dp. (We have dp = 190, 190, 495, and 924 for the 
nuclei ^"^Ne, ^^Ne, ^^Mg and ■^^Si, respectively). The sin- 
gular values decrease with increasing index and rapidly 
become exponentially small. (Degeneracies are due to 



spin/isospin symmetry). This suggests that a truncation 
of the sum in Eq. |(2Jl should yield an accurate approxi- 
mation to the ground-state. In fact, the density matrix 
renormalization group (DMRG) ||2J|, exploits this rapid 
fall-off of singular values in a wave-function factoriza- 
tion. For obvious reasons, the expansion (j^J is a factor- 
ization, and the correlated proton and neutron states are 
the factors. In what follows, we will devise a method 
that directly obtains the most important factors without 
knowledge of the exact ground-state. 

B. Derivation of main results 

To determine the optimal factor for a given truncation 
we make the ansatz 

IV') = Elft)K) (3) 

for the ground state. Here, the unknown factors are the 
proton-states \pj) and the neutron-states \nj). These 
states may be correlated and need not to be normalized. 
The truncation is controlled by the parameter fl which 
counts the number of desired factors. Figure ^ suggests 
that <C min((ip,djv) yields accurate approximations 
to shell- model ground states. This is also the result of 
our numerical computations below. 

Let H be the nuclear many-body Hamiltonian. To 
determine the unknown proton-states \pj) and neutron- 
states \nj) in Eq. ||2J) we consider the energy E = 

{■>p\H\4')/{'4'\'4')- Its variation 6E = yields {j = 
n 

J2{{nj\H\n,)-E{n,\n,)) \p,) = 0, 
n 

Y,{{P,\m}-E{p,\p,)) \n,) = 0. (4) 

i=l 

The solution of these nonlinear equations determines the 
optimal factors and the ground-state energy. Note that 
for fixed neutron (proton) states the first (second) set of 
these equations constitutes a generalized eigenvalue prob- 
lem for the proton (neutron) states. To fully understand 
the structure of the matrices involved we rewrite the first 
set of the Eq. Q as 



f {ni\H\ni) {ni\H\n2) 
{n2\H\ni) {n2\H\n2) 

V {nn\H\ni) 



{ni\H\nn) \ j \p,) \ 
{n2\H\nn) \p2) 

{nn\H\nn) ) \\Pn) ) 



E 



I {ni\ni)Ip {ni\n2)Ip 
{n2\ni)ip {n2\n2)ip 

\ {nn\ni)ip 



{ni\nn)Ip \ 
{n2\nn)ip 

{nn\nn)ip J 



( \P^) \ 

\P2) 

\ \Pn) j 



(5) 
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Here, \pj) = {pj,i,pj^2, ■ ■ ■ ,Pj,dp)'^ is a dp-dimensional 
column vector denotes the transpose) while Ip denotes 
the identity matrix in proton space. Thus, 



{n,\nj)ip = 



( {ni\nj) 
{ni\nj) 

\ 



\ 





(6) 



is a diagonal constant matrix of dimension dp x dp. The 
proton-space operators {ni\H\nj) stem from the nuclear 
structure Hamiltonian 



H = H]\[ + Hp + VpN, 



(7) 



with 
Hn 

Hp 

VpN 



n n,n' ,m.rn' 

^Spdldp+j ^ Vpp,q,qdlal,dqaq', (8) 



p,p',q,q' 



Here, we use indices p, q and m, n to refer to proton and 
neutron orbitals, respectively. The antisymmetric two- 
body matrix elements are denoted as Vijki- 

Thus, the proton-space Hamilton operator {ni\H\nj) 

is 



C. Treatment of symmetries 

Most modern shell-model codes use a basis of Slater de- 
terminants that preserve axial symmetry. It is straight 
forward to include this symmetry into the method pro- 
posed in this work. For such an "m-scheme" ground-state 
factorization we modify the ansatz ^ as 



M f2„ 



(10) 



m=-M k=l 



Here \pk{m)),k = l,...,fi„i (|nfe(-m)),fc = l,...,f2„) 
denote dp,™ dimensional proton states ( dN,-m dimen- 
sional neutron states) with angular momentum projec- 
tion Jz = m {Jz = —m), and the sum over m runs over 
all possible values of Jz- The ansatz ((Tn|l leads to a gen- 
eralized eigenvalue problem similar to Eq. with the 
only difference that permissible products of proton states 
and neutron states have zero angular momentum projec- 
tion: 

M n„ 



J2 E ("'^'("Oi^i^'^M) 

m=-M k=l ^ 

^E{nk'{m')\nk{m))^ \pk{-m)) = 0, (11) 



{ni\H\nj) = E Vpnp'n' (rtj |a^a«' a],dp/ 

+ {ni\HN\nj) + {nt\nj)Hp. (9) 

Note that the neutron-proton interaction Vppf results 
into a one-body proton operator while the neutron 
Hamiltonian Hn yields a constant. This concludes 
the detailed explanation of the first set of equations in 
Eq. Q. The second set has an identical structure, only 
the role of neutrons and protons is reversed. 



J 



J2 J2((P^'('^')\H\Pkim)) 

m=-AI k=l ^ 

-E{pk'{m')\pk{m)) I \nk{~m)) 



(12) 



These equations have to be fulfilled for all possible values 
of m' and k' . It is again useful to display the eigenvalue 
problem 11111 for the proton states in more detail 



/{n{-M)\H\n{^M)) ■■■ {n{-M)\H\n{M)) \ / |p(M)) 

V {niM)\H\n{-M)) ••• {n{M)\H\n{M)) J \ \p{-M)) 

I 



Here, we used the shorthands 



E 



{n{-M)\n{-M)) \p{M)) 



{n{M)\n{M)) \p{-M)) 



(13) 



the rectangular block matrices 



|p(m)) = (Ipi(to)), |_p2(m)), • • ■ , |k2^(to)))^, (14) 



{n{k)\H\n{l)) 



4 



/ (ni(fc)|ff|ni(Z)) ■•■ {ni{k)\H\nn,{l)) \ 

\ {nnAk)\H\n,il)) ••• {nnAk)\H\nn,{l)) ) 

and the overlap matrices 

{n{k)\n{k)) = 

/ (ni(fc)|ni(fc)}/p^_fc ••• (ni(fc)|noj^(fc))/p__fc \ 

\ {nn,{k)\ni{k))ip^^k ■■■ (no, (fc)|nn, (fc))/p,-fc / 

Note that the right hand side of Eq. H13() is a vector. Note 
also that {ni{k)\H\nj(l)) is a rectangular dp^k x c?p,i ma- 
trix similar to Eq. JHl, while Ip^k denotes the identity 
operator for the proton sub-space with angular momen- 
tum Jz = k. The eigenvalue problem (|13|l differs from the 
eigenvalue problem (0) due to the block diagonal overlap 
matrix on its right hand sight. The eigenvalue problem 
for the neutrons is identical to Eq. when the role of 
neutrons and protons is reversed. 

The number of factors used in the m-scheme factoriza- 
tion is given by the parameters f2„i- These parameters 
are input to the factorization. In the following, we use 

f^m = ^m{a) = max(l,Q!(ip^m), (15) 

and recall that dp,m is the dimension of the proton- 
subspace with angular momentum projection — m. 
For a — the most severe truncation is obtained and 
leads us to solve eigenvalue problems of the dimension 
dp and djv, respectively. Setting a = 1 leads to an 
eigenvalue problem with the same dimension as an ex- 
act diagonalization in m-scheme. Below we will see that 
the choice (|15ll of parameters yields rapidly converging 
results. However, there may be other parameterizations 
that are superior. Comparison of Eq. and Eq. © 
shows that the m-scheme factorization yields a lower- 
dimensional eigenvalue problem than the factorization 
© if the same number of factors is used. 

Note also that a j j-coupled scheme can be used. Let 
\pk{J,m)) {\nk{J,m))) be a proton (neutron) state with 
angular momentum quantum number J and angular mo- 
mentum projection ~ m. A ground-state of an even- 
even nucleus has J = and can be factored as 

,7 nj,„, 

IV) =^ ^ ^ |pfc(J,m))|nfc(J,-m)), (16) 

J rn— — J k—1 

where fij^m is the number of retained states with angular 
momentum quantum number J and angular momentum 
projection = m. We choose 

i^j,rn{a) = max (1, adp,j,„0, (17) 

where dp^j^m is the dimension of the proton-subspace 
with angular momentum J and projection m. This num- 
ber is independent of m for fixed J. Generalizations of 
the ansatz (I16|l to nonzero J are straight forward. 



Other symmetries like parity can also be used to fur- 
ther reduce the dimensionality of the eigenvalue prob- 
lem. Parity even states, for instance, are products of 
parity even proton states with parity even neutron states 
or products off parity odd proton states with parity odd 
neutron states. The ability to implement symmetries is 
particularly important as it widens the flexibility of the 
ground-state factorization. Consider for instance shell- 
model problems with proton and neutron spaces that 
differ considerably in size such that dp <^ djv. In such 
cases one might switch to a more symmetric factorization 
and factorize the ground-state as follows: one of the fac- 
tor spaces would consist of neutron states that are based 
on a subset of neutron single-particle orbitals, while the 
other factor states would be based on the remaining neu- 
tron orbitals and all proton orbitals. In such a scenario, 
the correct implementation of the Pauli principle requires 
care. 



D. Numerical solution and computational details 

We discuss the solution of the eigenvalue problems. 
For notational convenience we focus on the solution of 
the equations Q). The corresponding equations Hll|) 
and (|12() of the m-scheme factorization or the can be 
treated similarly. We solve the coupled set of nonlin- 
ear equations Q in an iterative procedure. We choose 
a random set of linearly independent initial neutron- 
states {\nj), j = 1, . . . , r2} and construct the Hamiltonian 
and overlap matrix presented in Eq. We then solve 
this generalized eigenvalue problem of dimension fldp for 
those proton-states {\pj),j = l,...,r2} that yield the 
lowest energy E. In practice, we use the sparse matrix 
solver Arpack "33^ for this task. The solution of the gen- 
eralized eigenvalue problem requires us to provide the LU 
factorization of the overlap matrix (i.e. the right hand 
side of Eq. jSJ). This factorization simplifies consider- 
ably since the overlap matrix is a direct product Ip(^Mis[ 
and only requires the LU-factorization of the x di- 
mensional overlap matrix Mn with elements {rii \nj). We 
then input the resulting proton-states to the second set 
of Eq. The solution of this ^Idpj dimensional problem 
yields improved neutron-states {\nj),j = 1, . . . , $7} and 
an improved (i.e. lowered) ground-state energy E. We 
iterate this procedure until the ground-state energy E is 
converged. 

For small values of H. we typically need about 20 it- 
erations to obtain a converged energy, and the number 
of iterations decreases with increasing number of kept 
states fl. For maximal value il — d^ , the first solu- 
tion of the proton-problem ((SJ already yields the exact 
ground-state. This can be seen as follows: Using Slater 
determinants {\nj) = = 1, • ■ • , i^} as input to the 

proton-problem (|3J) yields a matrix-problem that is iden- 
tical in structure to the full shell-model problem. 

Let us compare the effort of the ground-state factoriza- 
tion with an exact diagonalization. Both methods require 
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all Hamiltonian matrix elements {■Ka\{vf3\H\v^)\'Ks). The 
advantage of the factorization method is that the dimen- 
sionality of the eigenvalue problem is x max((ip,djv) 
with r2 ^ dp,dN while the dimension of the exact di- 
agonalization scales like dp x d^- Note that existing 
shell-model codes may easily be modified to include the 
ground-state factorization in order to obtain accurate ap- 
proximations or to compute useful starting vectors for a 
Lanczos iteration. 

It is possible to reduce the Eqs. Q to a standard 
eigenvalue problem. For this purpose we choose a ran- 
dom orthonormal set of initial neutron-states {\nj),j ~ 
1, . . . , fl} as input to the first eigenvalue problem. This 
reduces the "overlap" matrix {nj\ni) to a unit matrix, 
and we solve a standard eigenvalue problem to obtain the 
proton-states {\pj),j = 1, • ■ ■ , The resulting proton- 
states will not be orthogonal since they are components 
of only one solution vector of dimension Qdp. Their coef- 
ficient matrix C with elements Caj = {na\pj) may, how- 
ever, be factorizcd in a singular value decomposition as 
C = UDV^. Here D denotes a diagonal x 17 matrix 
while [/ is a (column) orthogonal dp x fl matrix, and V 
is a orthogonal 17 x f2 matrix. The transformed states 

dp 

\p'j) = Xl^^ji'^"^' 

Q = l 

n 

In'; 



3' 



i=l 



are orthonormal in the proton-space and in the neutron 
space, respectively, and they fulfill 



(18) 



The orthogonal proton-states {\p'j),j = 1, ... ,17} are 
then input to the second set in Eq. which poses a 
standard eigenvalue problem for the neutron-states. Note 
that the transformed neutron states Djj\nj) are useful 
starting vectors for the Lanczos iteration of the sparse 
matrix solver. The resulting neutron-states should again 
be orthonormalized by a singular value decomposition. 
These singular value decompositions are very inexpensive 
compared to the diagonalization. Similarly, we may cast 
the generalized m-scheme eigenvalue problem (|13f) into a 
standard eigenvalue problem by enforcing orthogonality 
(nfe(m)|n/(m)) oc Sk,i between neutron states with iden- 
tical angular momentum projection through the singular 
value decomposition. 

In our implementation of the m-scheme factorization, 
we compute the matrices of the proton space operators 
Hp and d^cLp' as well as the matrices of the neutron space 
operators and aj^a„'. These sparse matrices can be 
stored in fast memory. The sparse matrix on the left hand 
side of Eq. (0 is constructed from these matrices. This 
sparse matrix can be kept in fast memory for sufficiently 



small problems; for larger problems, its nonzero matrix 
elements along with the row-column information can be 
stored in large chunks of data on disk without a severe 
increase of computing times. 

Our implementation of the jj-coupled factorization is 
somewhat more tedious. Starting from the m-scheme 
basis states \pk{m)) and |nfe(m)), we create basis states 
|pfe(J, m)) and \nk{J,m)) with good angular momentum 
by numerical projection. The projection operator is |3l| 



P.J 



n J(J+1)_^-(J- + 1)' 



(19) 



where the product over j runs over all possible angular 
momenta, and J denotes the angular momentum opera- 
tor. The matrices of the proton Hamiltonian Hp and the 
neutron Hamiltonian H^ are then transformed to this 
basis and stored. The matrices of the operators a^Op' 
and dl^dn' are not sparse in the basis with good J, and 
are therefore kept in the m-scheme basis. We perform 
the appropriate basis transforms in the construction of 
the Hamiltonian matrices. 



III. NUMERICAL RESULTS 

A successful application of the ground state factoriza- 
tion would yield accurate approximations from calcula- 
tions involving rather small dimensional eigenvalue prob- 
lems. Clearly, the outcome depends on the model space 
and the interaction under consideration. In this section 
we apply the m-scheme factorization to realistic struc- 
ture calculations involving sd-shell nuclei, p/-shell nuclei 
and no-core shell-model computations for '*He. The re- 
sults are compared to exact diagonalizations. We expect 
the method to converge most rapidly in the case of weak 
proton-neutron correlations. Thus, the study of N = Z 
nuclei provides a challenging testing ground since T = 
proton-neutron correlations may be strong in these sys- 
tems. Throughout this section d denotes the dimension 
of the m-scheme eigenvalue problems (|11|) and (|12|) we 
actually solve, while dmax denotes the m-scheme dimen- 
sion required for an exact solution of the problem. 



A. sd-shell nuclei 

We apply the m-scheme factorization to the sc?-shell 
nuclei ^'^Mg, ^^A\, ^^Mg, and ^^Al and use the USD inter- 
action. While the factorization is particularly suited to 
compute accurate approximations of the ground states, 
we may also use it for the computation of low-lying ex- 
cited states. In some cases, excited states can be ob- 
tained as a by-product of the ground-state computation. 
While solving the eigenvalue problems (|ll|l and (|12|l for 
the ground-state, we may also compute excited state so- 
lutions. Let us assume that we solve the equations (|llf) 
for the proton states. The excited proton state solutions 
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FIG. 2: Low-energy spectrum of sd-shell nuclei versus the 
dimension d of the eigenvalue problem. Results are obtained 
from targeting the ground state. Left panel: Mg isotopes 
(dmax = 28503) ^**Mg (top), ^"Mg (bottom). Right panel: Al 
isotopes (dmax = 69784) ^^Al (top) ^**A1 (bottom). 



will be obtained in the presence of neutron states that are 
optimized for the ground state. Therefore, we expect that 
the excited states are less accurately reproduced than the 
ground-state. Figure [21 shows the resulting low-energy 
spectra for "^^Mg, ^^Mg, ^^Al, and ^SAl. 

The ground states converge most rapidly as more fac- 
tors are retained in the factorization and the dimension d 
of the eigenvalue problem increases. Typically, excellent 
results are obtained from computations involving relative 
dimensions d/dmax ~ 1/4... 1/3. For the Mg isotopes, 
excited states converge somewhat slower than the ground 
states, and level spacings are reproduced to a very good 
accuracy already at severe truncations. This shows that 
the factors of the low- lying excitations have a large over- 
lap with the corresponding ground state factors, and we 
assume that this finding is related to the band structure 
of the low-lying excitations. The situation is different 
for the Al isotopes, as evident from the right panel of 
Fig. 121 This slower convergence of excited states is not 
unexpected due to the absence of band structure in these 
nuclei. 

There are at least two approaches to improve the con- 
vergence of the excited states. In the first approach, we 
may target excited states by solving the eigenvalue prob- 
lems (|llll and (|12|l directly for an excited state. This 
method works well for the first and second excited states 
of ^"'Mg, as shown in Fig. |2| However, this approach is 
unstable for higher excited states of ^^Mg and for the first 
excited state of ^^Al, as the solutions of the eigenvalue 
problem are oscillating but fail to converge with increas- 
ing number of iterations. In the case of ^^Al we proceed 
as follows. The ground state has angular momentum 
J = 5 while the first excited state has angular momen- 
tum J = 0. We may thus use the jj-coupled ansatz Hlt)|) 




30000 



FIG. 3: . Energies of the first and second excited states of 
^*Mg versus the dimension of the eigenvalue problem. Hollow 
data points: results from targeting the ground state. Filled 
data points: results from directly targeting the excited states. 
Dashed lines: exact results. 



to compute the lowest-lying state with angular momen- 
tum J = 0. This yields the first excited state. Figure ^ 
shows the results plotted versus the corresponding m- 
scheme dimension. The convergence is much improved 
and comparable to that of the ground-state. 

So far we have only focused on the energies. In the 
remainder of this section we discuss how the states and 
their quantum numbers are reproduced by the factoriza- 
tion. For ^"^Mg we analyze the wave-function structure 
of the low- lying levels. Figure [S] shows the squared over- 
laps with the exact results, obtained from targeting the 
ground state. The solution of an eigenvalue problem of 
only 10% of the full dimension dmax already yields be- 
tween 90-96% of squared overlap. The directly targeted 
ground-state is reproduced to more than 99% once the 
dimension exceeds 20% of the full dimension dmax- The 



inset of Fig. El shows that the defect 1 



(V'exactlV'factor)^ 
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FIG. 4: Energy of the first excited state of ^''Al versus the 
dimension of the eigenvalue problem. Hollow data points: 
results from targeting the ground state. Filled data points: 
ground state calculation for zero angular momentum. Dashed 
hue: exact result. 
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FIG. 5: Squared overlaps (t/'oxact I'/'factor)^ for the low-lying 
states in ^'*Mg versus the dimension of the eigenvalue problem 
(from targeting the ground state). The inset shows the devia- 
tion 1 — (V'oxact I '(/'factor)'^ versus the dimension of the eigenvalue 
problem. 



FIG. 6: Angular momentum value j for the low-lying states 
in ^*Mg versus the dimension of the eigenvalue problem (from 
targeting the ground state). 



B. p/-shell nuclei 



decreases exponentially fast as more factors are retained. 
While the convergence is best for the directly targeted 
ground-state, the excited states are also very accurately 
reproduced. 

Let us also verify for ^''Mg that the quantum numbers 
of the low-lying states are reproduced correctly. Due to 
rotational symmetry, the expectation value for the angu- 
lar momentum should fulfill 

(i')=j(j + i), (20) 

where j is a nonnegative integer. Figure shows the j- 
values of the low-lying states for ^''Mg. The results were 
obtained by targeting the ground state. The angular- 
momenta are very accurately reproduced once about 20% 
of the states are retained in the factorization. This is re- 
markable since rotational symmetry is not enforced in 
the m-scheme factorization, and no kind of constraint 
or projection was used. Similar results are obtained for 
the total isospin and for other sd-shell nuclei. The accu- 
rate reproduction of quantum numbers, wave-functions 
and energies implies that transition matrix elements can 
accurately be computed. This concludes or detailed dis- 
cussion of ^"'Mg. 

We finally mention that we also compared the m- 
scheme factorization with the jj-coupled scheme. For 
^^Mg we find practically identical convergence of the 
ground-state energy when plotted versus the relative di- 
mension d/dniax of the corresponding eigenvalue problem. 
The dimensions d, c?max of the eigenvalue problem in the 
jj-coupled scheme are, of course, smaller than for the 
771-scheme. However, this does not translate directly into 
a computational speed-up since the jj-scheme algorithm 
is more complex and involves less sparse matrices. We 
believe that its main advantage consists of the possibility 
to directly target the lowest state with a given angular 
momentum. 



Many theoretical results for p/-shell nuclei are avail- 
able from exact calculations for the KB3 interaction. In 
the lower p/-shell, diagonalizations can be based on an 
TO-scheme basis 0]. The m-scheme dimensions of upper 
p/-shell nuclei exceed one billion, and exact diagonaliza- 
tions have been performed in a J = coupled basis 3|, 
reducing the dimensions to the order of ten million. In 
this section we compare the results from m-scheme fac- 
torization with the exact results. We are particularly 
interested in the efficiency of the method and would like 
to answer the following question. How does the relative 
dimension cJ/dmax, at which an accurate approximation 
to the ground state is obtained, scale with increasing di- 
mension c?inax required for an exact diagonalization? 

For p/-shell nuclei we use the KB3 interaction [3^ . 
Figure[7|shows the low- lying energies for *^Cr plotted ver- 




FIG. 7: Data points; Low-lying energies of ''^Cr (KB3 in- 
teraction) versus the dimension d of the eigenvalue problem 
relative to the m-scheme dimension dmax. Dashed lines: Ex- 
act results. Dashed-dotted line: Exponential fit _E(d/dniax) = 
—32.92 + 0.851 exp (—31.38 d/dmax) to the ground state en- 
ergy. 
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FIG. 8: Angular momentum value j for the low-lying excita- 
tions of the p/-shell nucleus ''^Cr (KB3 interaction) plotted 
versus the dimension d of the eigenvalue problem relative to 
the m-scheme dimension dmax. The dotted lines are the exact 
results. 



> 

. \ Si 

\ ^ 
\ ^ 

\ +0.1 


1 1 1 1 1 1 ■ 1 1 1 


- 


Gl ^^-^^Ve results 

-V -- 0.576 exp(-236 d/ d^j^j,) - 






0.003 0.006 0.009 
^ d/d 

max 




— °''^Ni results 
■ — - 0.337 exp(- 1 392 d / d^^^^) 


1 .... 1 . 



0.0005 0.001 



d/d 

max 

FIG. 9: Ground-state energy E versus the dimension d of 
the eigenvalue problem relative to the m-scheme dimension 
dmax ^ 1.09 X 10^ for ^"^Ni. The data points are from the m- 
scheme factorization, and the dashed line is an exponential fit 
to the data. Inset: Similar plot for ^"Fe (dmax « 110 x lO'*). 



sus the relative dimension d/dma.^ of the eigenvalue prob- 
lem. The exact ground state energy is Eq = —32.95 MeV 
and results from the solution of an eigenvalue problem 
with dimension dmax — 1-96 x 10^ 0. The ground-state 
energy converges exponentially quickly as the number of 
retained factors increases. The right-most data point re- 
sults from an eigenvalue problem with only 8% of the 
m-scheme dimension and involves fl = ^^^^ f2„, — 391 
factors. It deviates less than 100 keV from the result 
of an exact diagonalization. An exponential fit of the 
form E{d/dn^iix) = a-|-&exp {^cd/ d-a^a-x) to the right-most 
six data points yields the estimate E{\) = —32.92 MeV, 
which is only 30 keV above the exact result. The excited 
states are obtained from targeting the ground state. The 
level spacings of the two lowest excitations are very ac- 
curately reproduced even at the most severe truncation, 
while the spacings to the higher levels are about 200-300 
keV too large. The angular momentum expectation val- 
ues are plotted in Fig.|Sl For d/dmax 0.04, energies and 
quantum numbers are sufficiently well converged. Con- 
sidering the modest size of the eigenvalue problem we 
solved, these are very good results. 

After this detailed discussion of ^*Cr, we factor the 
ground states of ^°Fe and ^^Ni. The exact energies 
are E = -67.0MeV and E = -78.46MeV, respec- 
tively 01, and the corresponding m-scheme dimensions 
are dmax ~ HO x 10^ and dmax ~ 1-09 x 10^, respectively. 
Fig. 1^ shows that the ground states of these nuclei can 
very efficiently be factored. Using an exponential fit to 
the numerical data points, the ground state energies are 
reproduced within a deviation of 50 keV for ^°Fe and 
100 keV for ^^Ni. Most importantly, the relative dimen- 
sion of the eigenvalue problem we solve is d/dmax ~ 1% 
for ^°Fe (from = = 352 states) and about 

d/dmax « 0.1% for (from n = = 147 

states) . This suggests that the factorization gets increas- 
ingly efficient as the dimension of the problem increases. 



C. No-core shell model 

As a final test, we consider a no-core shell- model prob- 
lem and apply the m-scheme factorization to ^He using a 
manageable model space and a realistic interaction from 
a G-matrix calculation. The model space consists of the 
0s-0p-ls-0d-0/-lp shells. The G-matrix stems from a 
Ibhio calculation and is based on the Idaho-A potential 
[3^ . The Idaho-A potential is derived from an effective 
Lagrangian that respects QCD inspired chiral symmetry. 
We calculate the G-matrix from 

GH = V+ ?^G(c^) (21) 

uj - QTQ 

where uj is the starting energy, T is the kinetic energy 
operator, and Q is the Pauli operator. Our Pauli oper- 
ator allows for all allowed configurations to be active in 
the chosen model space. We also employ folded-diagrams 
calculated at a) = —20.0 MeV to decrease the dependence 
of the resulting two-body interaction on the starting en- 
ergy. Details concernin g th e derivation of the G-matrix 
may be found in Ref. |33- We also include center-of- 
mass corrections perturbatively. Finally, we employ the 
method of Ref. ^33| to obtain an interaction that yields a 
ground-state energy that is approximately free of center- 
of-mass contamination. Note that this small space is not 
sufficient to completely describe the ^He nucleus, but still 
illustrates the power of the factorization method for prob- 
lems in which the core is absent. 

Figure [TUI shows the energies for the three lowest lying 
states of *He versus the dimension d of the eigenvalue 
problem we solve. The results for the excited states were 
obtained while targeting the ground state. The exact re- 
sults are obtained from a diagonalization with m-scheme 
dimension dmax — 79298. Note that the ground state 
and the excited states converge very fast toward the ex- 
act results. A calculation with d/dmax ~ 0.2 already 
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yields excellent approximations, and the angular momen- 
tum quantum numbers are converged. 




80000 



FIG. 10: Energies of low-lying states of ^He plotted versus 
the dimension d of the eigenvalue problem. The dashed lines 
are the exact results. 

Let us finally suggest an alternative treatment of the 
center of mass problem. The center of mass is separable 
in an oscillator basis where all many-body states with up 
to A^max oscillator quanta are included. The factoriza- 
tion could be applied in this scheme by combining proton 
states with n oscillator quanta and neutron states with 
n' oscillator quanta such that n + n' < A^max- 



IV. CONVERGENCE OF THE METHOD 

A. Convergence properties 

The results of the preceding section showed that the 
factorization converges exponentially quickly as more fac- 
tors are retained. So far we considered nuclei with equal 
dimension of proton space and neutron space, most of 
them being N = Z nuclei. What can be expected for 
other cases? To answer this important question, we com- 
puted the exact ground states of several p/-shell nuclei 
and numerically performed singular value decompositions 
of their amplitude matrices '^ap as defined in Eq. 
Fig. ^2 shows logarithmic plots of the resulting singu- 
lar value spectra. The singular value spectra exhibit a 
very sharp initial falloff followed by an exponential de- 
cay. The initial falloff is stronger for larger dimension 
of the proton space dp, and this renders the factoriza- 
tion method very effective. There is no clear trend for 
isotopic chains. The singular value spectra of the lighter 
p/-shell nuclei decay most rapidly for the N = Z nuclei, 
while the decay is faster for mid-shell nuclei away from 
N = Z. Our results suggest that the application of the 
factorization method is not limited to N = Z nuclei. We 
also computed the number of factors, ^{x), such that 
^^i^) „2 ^ ^ f_ ^ ^ 0.99 and x = 0.999. The results 




FIG. 11: Singular value spectra (squared singular values s|) 
for ground states of Ti isotopes (upper left, dp — 190), V 
isotopes (upper right, dp — 1140), Cr isotopes (lower left, 
dp = 4845), and ^"Mn (lower right, dp = 15504). 



of Table ^ show that the factorization of all these nu- 
clei should converge rapidly as more factors are included. 
Larger dimensional model spaces usually require the re- 
tention of more factors, but the ratio fl(x)/ min (dp, (In) 
decreases with increasing size of the problem. Note that 
^*Cr is relatively difficult to factor into proton and neu- 
tron states. This suggests that this nucleus exhibits par- 
ticularly strong proton- neutron correlations. 

We recall that the m-scheme factorization (|1U|I requires 
as input the number of factors with a given angular mo- 
mentum projection, Q„n which were taken according to 
Eq. H15|) . It is interesting and important to check this 
choice of input parameters. To this purpose we com- 
pare the singular value spectrum from the factorization 
with the singular value spectrum from an exact calcu- 
lation. The factorization was performed for ^*Cr using 
a = 0.04 and a — 0.08 in Eq. ifT^ . These truncations 
included a total of f7 = Y.m ^(") = 19"^ and f2 = 391 



Nucleus 


dp 


djv 


J7(0.99) 


J7(0.999) 




190 


190 


40 


75 


45 


190 


1140 


45 


109 


«Ti 


190 


4845 


43 


126 


46y 


1140 


1140 


98 


270 


47y 


1140 


4845 


111 


324 


48y 


1140 


15504 


151 


392 


4«Cr 


4845 


4845 


258 


775 


48 Cr 


4845 


15504 


168 


619 


50Mn 


15504 


15504 


197 


821 




15504 


15504 


163 


570 



Sj>x for 



TABLE I: 

Proton space dimension dp and neutron space dimension djv 

for various p/-shell nuclei. Q{x) denotes the number of 
factors that have to be retained for an overlap X/j^-i •^j ~ ^ 
with the exact ground state. 



10 



factors, respectively. Figure El compares the resulting 
singular value spectra with the singular value decompo- 
sition of the exact ground state. The agreement between 
the exact results and the approximations is rather good, 
and improves with increasing number of retained factors. 
Note however, that the smaller singular values deviate 
from each other. This suggests, that it should be possi- 
ble to somewhat improve the choice of input parameters 

We finally mention that we lack an understanding 
about the rapid initial falloff in singular value spec- 
tra. For the Hamiltonians considered in this work the 
falloff is evident. Moreover, the results obtained from 
DMRG calculations over the past decade demonstrate 
that density-matrix spectra decay rapidly for the ground 
states of a large number of relevant Hamiltonians. A few 
works address the theoretical foundations of the DMRG. 
Peschel and coworkers investigated density matrix spec- 
tra of several soluble problems j3Zj, while Okunishi et 
al. discussed the asymptotic behavior of density matrix 
eigenvalues for noncritical spin systems [33 |. Ostlund 
and Rommer showed that if the DMRG renormaliza- 
tion converges to a fixed point, the DMRG ground-state 
is of a special matrix-product form [s^- Given the 
lack of generally valid analytical results, it is thus in- 
teresting to numerically investigate singular value spec- 
tra for "generic" Hamiltonians. For this purpose we 
considered the model space of p/-shell nuclei ^''Ti and 
"^^V and used a random two-body interaction that pre- 
serves spin and isospin, i.e. the spin/isospin coupled 
two-body matrix elements Va.p are independent Gaus- 
sian random variables with zero mean {Va^p) = and 
variance {Va^pVa-^p-) = {5a,a'5p^p- +5a,p'5p,a')- For a re- 
alistic choice of the single-particle energies we find singu- 
lar value spectra that are similar to the realistic spectra. 
However, setting all single-particle energies to zero yields 
singular value spectra with longer tails and a less rapid 
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index j 

FIG. 12: Singular value spectra (squared singular values s|) 
for ■'^Cr. Singular value decomposition of exact ground state 
(full line), m-scheme factorization using Q. = 391 factors 
(dashed line), and n = 197 factors (dashed-dotted line). 



decay. 



B. Comparison with other truncation methods 

In recent years, several truncation methods have been 
developed for and applied to shell model problems. In 
this subsection we compare some of these approaches 
with the method presented in this work. This comparison 
focuses on convergence and accuracy of low-lying energy 
spectra at a given level of truncation. 

We start with the DMRG which also bases its trunca- 
tion on the singular values |24j . Dukelsky and coworkers 
applied the DMRG to nuclear structure problems involv- 
ing pairing |26l | and pairing-plus-quadrupole interactions 
|27|. These applications were very successful as accu- 
rate results could be obtained for huge Hilbert spaces. 
The factorization method proposed in this work can only 
treat much smaller Hilbert spaces and cannot compete 
with the DMRG for these systems. However, the recent 
DMRG calculation "28] for the realistic nuclear structure 
problem (^"^Mg with USD interaction) converges very 
slowly. 

Recently, Andreozzi et al. T^,'20'| used a small number 
of correlated proton states and correlated neutron states 
as a truncated basis for shell model problems. The cor- 
related proton (neutron) states are the low-energy eigen- 
states of the proton-proton (neutron-neutron) Hamilto- 
nian. The full shell model Hamiltonian including the 
proton-neutron interaction is then solved in this space. 
Andreozzi and Porrino report exponentially converging 
results and a considerable reduction in the number of ba- 
sis states p^ . This procedure differs from our approach 
mainly by the absence of a variational principle. 

A third related method is the exponential conver- 
gence method (ECM) developed by Horoi and coworkers 
111 III III 

llq . In this method, shell- model configu- 
rations are ordered according to their average centroid, 
which are obtained from statistical spectroscopy. This or- 
dering gives a natural truncation scheme, and analytical 
arguments suggest an exponential convergence of ener- 
gies with increasing number of retained configurations. 
Once the exponential region is identified, the full space 
energies can be extrapolated by an exponential fit. A di- 
rect comparison is not easy since the FPD6 interaction is 
used for p/-shell nuclei, and since ECM results are plot- 
ted versus JT-coupled dimension of the truncated space. 
We believe, however, that our method converges more 
rapidly than the ECM. For ^^Cr, for instance, our rate of 
exponential convergence is c « —31.38 (See Fig.O, which 
is about a factor eight larger than what is reported for 
the ECM in Fig. 1 of Ref. |l||. For ^'^Ni, our exponential 
rate is about a factor 200 larger than the ECM rate 4(| , 
and our identification of the exponential region requires 
a m-scheme dimension d ~ 10® (See Fig. I10() while the 
ECM requires an m-scheme dimension of 4-5 million 40J . 

For upper p/-shell nuclei, truncations can be based 
on the maximal number t of nucleons outside the /7/2 



11 



subshell [4l|. Within this truncation, the convergence of 
the energy is rather slow, and it is difficult to extrapo- 
late from results in truncated spaces to the full Hilbert 
space. Mizusaki and Imada devised extrapolation meth- 
ods that link the error due to the truncation to the vari- 
ance of the energy in the truncated space. This approach 
lead to a first order ^V?\ and a second-order extrapola- 
tion method 18] for predictions of low-lying states in 
various p/-shell nuclei. For ^^Ni the approximation of 
a closed fj/2 subshell is well justified, and the extrap- 
olation methods yields results that are superior to the 
factorization [Z^]. For ^®Cr, however, the factorization 
seems to be of advantage: The exact ground-state energy 
being E = —32.95 MeV and the m-scheme dimension 
d = 1.96 X 10^. The first-order extrapolation method 
[13 yields E = -33.008 for t = 5 and E = -32.975 
for t = 6, and the corresponding m-scheme dimensions 
are d{t = 5) « 1.3 x 10^ and d{t = 6) « 1.76 x 10^. 
The second-order extrapolation ("scheme I") yields 
E — —32.91 for t — 5. Slightly better results are ob- 
tained from a different truncation scheme ( "scheme 11" ) . 
The ground-state factorization yields a comparably good 
energy estimate E = -32.92 MeV (See Fig. [TJ from 
solving a much smaller eigenvalue problem of dimension 
d = 1.6 X IQS. 

The mixed-mode shell model approach developed by 
Gueorguiev et al ,21^ combines single-particle config- 
uration and iS'C/(3) collective configurations to describe 
the interplay and competition between single-particle and 
collective degrees of freedom. For the sd-shell nucleus 
^■^Mg, the mixed-mode shell model yields good approx- 
imations to the binding energy (within 2% deviation of 
the exact result), and low-energy configurations which 
exceed 90% overlap with the exact results. These results 
stem from a truncated space of only 10% the full Hilbert 
space 1121 ■ At the 10% level of the truncation, the fac- 
torization method yields an energy deviation of less than 
1% (See Fig.[2J), and squared overlaps exceed 96% for the 
two lowest lying states and 90% for the following three 
states (See Fig. 01. 

The method proposed in this work thus compares well 
to most of the alternatives regarding convergence and ac- 
curacy at a given level of truncation. Note, however, that 
its implementation seems somewhat more complex than 



a shell-model approach with a configuration truncation 
and somewhat less complex than the DMRG algorithm. 



V. CONCLUSION 

We approximated the ground states of realistic nuclear 
structure Hamiltonians by sums over products of corre- 
lated proton states and correlated neutron states. The 
optimal states are determined by a variational principle 
and are the solution of rather low-dimensional eigenvalue 
problems. Computations for sd-shell nuclei, p/-shell nu- 
clei, and no-core shell models show that the method con- 
verges exponentially quickly as more factors are included, 
and that accurate approximations to shell-model ground 
states and low-lying excitations may be obtained. For 
the largest problems we considered, the dimension of the 
eigenvalue problem was reduced by three orders of mag- 
nitude. Momentarily, the application of this method is 
limited by the size of the proton space and the neutron 
space. An interesting future development would also con- 
sider the factorization of these spaces in order to treat 
larger dimensional problems. While the reason of the ex- 
ponential convergence is not yet understood, computa- 
tions of shell model problems with realistic and random 
two-body interactions suggest that this behavior can be 
expected for a variety of interactions. 
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